#include "AerosolConcentration.h"
#include "AtmosModel.h"
#include "common.h"

/**********************************************************************c
c              aerosol model (concentration)                           c
c              ----------------------------                            c
c                                                                      c
c                                                                      c
c  v             if you have an estimate of the meteorological         c
c                parameter: the visibility v, enter directly the       c
c                value of v in km (the aerosol optical depth will      c
c                be computed from a standard aerosol profile)          c

c                                                                      c
c  v=0, taer55   if you have an estimate of aerosol optical depth ,    c
c                enter v=0 for the visibility and enter the aerosol    c
c                optical depth at 550                                  c
c                                                                      c
c  v=-1          warning:  if iaer=0, enter v=-1                       c
c                                                                      c
c**********************************************************************/
void AerosolConcentration::parse(const long int _iaer, const AtmosModel& atms)
{
    iaer = _iaer;

	taer55 = 0.f;
	cin >> v;
	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */

	if(v == 0)
	{
		cin >> taer55;
		cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
		v = (float)(exp(-log(taer55/2.7628f)/0.79902f));
	}
	else if(v > 0) oda550(v, atms);
}

void AerosolConcentration::oda550(const float v, const AtmosModel& atms)
{
	/* aerosol optical depth at wl=550nm */
	/* vertical repartition of aerosol density for v=23km */
	/*                ( in nbr of part/cm3 ) */
	static const float an23[34] = {
		2.828e+03,1.244e+03,5.371e+02,2.256e+02,1.192e+02,
		8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
		5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
		4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
		2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
		7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
		1.078e-02,5.550e-05,1.969e-08,0.000e+00
	};


	/* vertical repartition of aerosol density for v=5km */
	/*                ( in nbr of part/cm3 ) */
	static const float an5[34] = {
		1.378e+04,5.030e+03,1.844e+03,6.731e+02,2.453e+02,
		8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
		5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
		4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
		2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
		7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
		1.078e-02,5.550e-05,1.969e-08,0.000e+00
	};

	taer55 = 0;
	if(fabs(v) <= 0) return;
	if(iaer == 0) return;

	for(int k = 0; k < 32; k++)
	{
		float dz = atms.z[k+1] - atms.z[k];
		float az = (115.f / 18.f) * (an5[k] - an23[k]);
		float az1 = (115.f / 18.f) * (an5[k+1] - an23[k+1]);

		float bz = (5.f * an5[k] / 18.f) - (23.f * an23[k] / 18.f);
		float bz1 = (5.f * an5[k+1] / 18.f) - (23.f * an23[k+1] / 18.f);

		float bnz = az / v - bz;
		float bnz1 = az1 / v - bz1;

		float ev = (float)(dz * exp((log(bnz) + log(bnz1)) / 2));
		taer55 += ev * sigma * 1.0e-03f;
	}
}

void AerosolConcentration::print()
{
	/* --- aerosol model (concentration) ---- */
	Output::Begin();
	Output::End();
	if(iaer == 0) return;

	Output::Begin();
	Output::Repeat(10, ' ');
	Output::Print(" optical condition identity :");
	Output::End();
	if(fabs(v) <= xacc)
	{
		Output::Begin();
		Output::Repeat(15, ' ');
		Output::Print(" user def. opt. thick. at 550nm :");
		ostringstream s;
		s.setf(ios::fixed, ios::floatfield);
		s << setprecision(4);
		s << setw(11) << taer55 << ends;
		Output::Print(s.str());
		Output::End();
	}
	else
	{
		Output::Begin();
		Output::Repeat(15, ' ');
		Output::Print(" visibility :");
		ostringstream s1;
		s1.setf(ios::fixed, ios::floatfield);
		s1 << setprecision(2);
		s1 << setw(8) << v << ends;
		Output::Print(s1.str());
		Output::Print(" km  opt. thick. 550nm :");
		ostringstream s2;
		s2.setf(ios::fixed, ios::floatfield);
		s2 << setprecision(4);
		s2 << setw(9) << taer55 << ends;
		Output::Print(s2.str());
		Output::End();
	}

	Output::Begin();
	Output::End();    
}


AerosolConcentration AerosolConcentration::Parse(const long int iaer, const AtmosModel& atms)
{
	AerosolConcentration aerocon;
	aerocon.parse(iaer, atms);
	return aerocon;
}

